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I.  STATEMENT  OF  THE  PROBLEM 


The  goal  of  this  project  is  to  improve  the  accuracy  of  rate  control 
systems.  The  traditional  direct-fire  gun-type  turret  has  an  accuracy 
requirement  on  the  order  of  0.2  milliradians  which  is  sufficient  for 
ranges  of  from  one  to  two  kilometers.  Laser  designators  and  guided 
missile  directors  are  designed  for  twice  the  range  and  therefore  require 
twice  the  accuracy,  i.e.,  approximately  0.1  milliradians.  Even  the  best 
modern  turret  systems  have  trouble  achieving  this  order  of  accuracy  all 
the  time.  This  effort  attempts  to  determine  which  system  character¬ 
istics  are  limiting  tracking  accuracy.  It  will  then  be  possible  to 
write  specifications  in  terms  of  hardware  characteristics  rather  than 
in  terms  of  system  performance  goals. 

The  premise  of  this  approach  is  that  the  solution  can  be  found  in 
the  nonlinear  behavior  of  the  turret  systems  at  very  low  rates.  The 
tracking  rates  of  interest  and  changes  in  rate  that  are  called  for  to 
achieve  precise  tracking  are  very  small.  A  30-kilometer-per-hour  target, 
for  example,  at  a  range  of  A  kilometers  and  a  heading  of  30  degrees 
would  result  in  a  1-milliradian-per-second  crossing  rate.  The  gunner 
will  command  small  changes  about  this  nominal  rate  in  an  attempt  to 
reduce  the  tracking  error.  We  can  infer  the  magnitude  of  these  changes 
to  be  less  than  0.3  milliradians  per  second  (referred  to  the  output)  by 
the  following  logic.  Since  the  desired  error  tolerance  is  on  the  order 
of  0.1  milliradians  the  gunner  should  have  the  ability  to  make  commands 
at  his  bandwidth  (3  radians  per  second)  which  would  result  in  amplitudes 
of  0.1  milliradian.  Assuming  a  sinusoidal  input  of  0.3-milliradian-per- 
second  amplitude  and  3-radian- per- second  frequency  the  output  would  be  a 
sinusoid  with  0 . 1-milliradian  amplitude. 

The  nonlinear  elements  in  the  turret,  i.e.,  backlash,  coulomb 
friction,  and  deadspace,  will  have  more  effect  on  the  turret  response 
at  low  rates  and  for  small  commands.  The  gunner  will  notice  a  decrease 
in  gain  and  an  increase  in  phase  lag.  He  will  attempt  to  compensate 
for  the  turret  changes  within  limits,  but  eventually  the  tracking 
error  will  increase  above  what  it  would  be  for  a  linear  system. 
Moreover,  the  human's  ability  to  compensate  will  vary  considerably 
with  training  and  between  people.  The  result  will  be  unpredictable 
system  performance  due  to  operator  differences  rather  than  hardware 
differences . 

Traditionally  the'  turret  spec  if icat ions  for  low  rate  tracking 
have  been  system  performance  specifications  with  a  human  operator  in 
the  loop,  e.g.,  0. 1-milliradian  root-mean- square  error  when  tracking 
a  1 .O-milliradian-per-second  target  for  10  seconds.  It  would  probablv 
be  better  to  specify  the  accuracy  requirements  against  sinusoidal 
inputs  and  better  still  to  specify  the  system  characteristics  without 
a  human  operator  in  the  loop.  Enough  analysis  of  linear  control 
systems  with  human  operators  has  been  done  to  relate  accuracy  to 
system  characteristics  (gain  and  phase  lag)  and  input  power  spectra  at 
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least  for  linear  control  systems.  This  report  addresses  the  system 
performance  at  low  rates  when  the  nonlinear  system  characteristics  are 
important.  Ultimately  this  work  should  lead  to  specifications  of 
system  characteristics  which  are  required  to  achieve  any  specified 
degree  of  accuracy. 


II.  APPROACH 

The  overall  solution  to  the  problem  requires  three  distinct  tasks 
to  be  carried  out  by  three  different  organizations. 

a.  First,  there  is  the  characterization  of  a  typical  turret 
control  system  at  very  low  rates.  An  M60A3  turret  control  is  being 
used  because  it  is  available  now  and  in  the  future  and  it  is  typical 
of  current  military  turret  control  systems.  Frequency  response  and 
transient  response  will  be  measured  at  various  amplitudes  from  low 
rates  near  or  below  threshold  up  to  high  rates  in  the  linear  region. 
This  work  is  being  done  under  contract  by  General  Electric, 
Pittsfield,  Massachusetts. 

b.  Second,  there  are  experiments  with  humans  in  the  loop  and 

with  computer  simulated  turret  response.  These  tests  will  charac¬ 
terize  the  human  operator  and  they  will  allow  for  some  limited  para¬ 
meter  variations  of  turret  response.  Human  Engineering  Laboratory  is 
doing  this  work.  4 

c.  Third,  there  is  an  analysis  task  with  both  the  turret 
response  and  the  human  operator  simulated  by  a  computer.  This  task 
allows  for  parametric  variations  of  turret  response  just  as  with  the 
human  experiments;  but  the  computer  allows  for  almost  unlimited  para¬ 
meter  variations  because  the  trials  are  faster  and  they  are  without 
random  human  variation.  This  report  describes  the  Phase  I  effort  on 
the  third  task. 

The  Phase  I  effort  develops  the  methods  and  computer  codes  based 
on  assumed  system  characteristics.  The  Phase  II  effort  will  do  it  all 
again  but  with  accurate  system  data  from  the  other  tasks. 

The  main  task  of  Phase  T  was  to  develop  a  self-optimizing  human 
operator  model.  The  model  must  minimize  root  mean  square  (rms) 
tracking  error,  subject  to  constraints  on  human  behavior,  in  a 
consistent  and  rational  manner.  This  model  was  then  used  to  determine 
the  effects  of  various  kinds  of  nonlinearities  on  system  performance. 
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III.  PROCEDURE  FOR  MAN-MODEL  OPTIMIZATION 


The  Phase  I  effort  had  a  goal  of  developing  techniques  for 
evaluating  nonlinearities.  This  included  the  man  model,  the  adaptive 
algorithm  for  the  man  model,  models  of  nonlinearities  and  some  limited 
data  to  show  how  well  the  models  work.  Complicated  turret  dynamics 
were  not  important  at  this  stage;  consequently  a  very  simplified  model 
of  turret  response  was  incorporated.  Figure  1  shows  a  block  diagram 
of  the  control  loop  and  the  linear  models  that  were  used  for  the  man 
and  turret.  The  reference  signal  was  either  white  noise,  a  sinusoid, 
or  a  maneuvering  tank.  The  human  model  was  a  conventional  linear 

model  from  Sheridan^  to  which  a  noise  remnant  was  added  for  reasons 
that  will  be  explained  later.  The  turret  response  was  given  a  time 

constant  of  0.1  seconds  and  a  gain  of  0.01  radians  per  second  per 
radian.  These  parameters  are  about  right  for  a  tank  turret  at  very 

low  rates.  Nonlinearities  were  introduced  in  the  digital  simulation 
at  the  points  shown  in  Figure  1. 

The  notation  in  Figure  1  was  chosen  to  be  consistent  with  the 

digital  simulation  shown  in  its  entirety  in  the  Appendix.  The  linear 
transfer  functions  were  simulated  by  the  step  invariant  zeta  transform- 
method. 


Figure  1.  Block  Diagram  of  Control  Loop  . 


^Man-Machine  Systems,  Sheridan  and  Ferrell,  1974,  MIT  Press. 
^Digital  Signal  Analysis,  Stearns,  1Q75,  Hayden  Book  Company  Tpc. 


An  algorithm  for  optimizing  the  human  transfer  function  was 
developed  while  using  white  noise  for  a  reference  signal.  It  minimized 
rms  tracking  error  but  with  a  penalty  for  open  loop  gain  margin  less 
than  6  db  and  open  loop  phase  margin  less  than  45  degrees.  Secant 
functions  accomplished  the  penalty  by  giving  a  cost  factor  of  one  at  6 
db  and  45  degrees  and  a  cost  factor  of  infinity  at  0  db  and  0  degrees. 

secant  {15  (Gain  Margin  -6)],  secant  (2  (Phase  Margin  -45)1 

The  loop  was  warmed  up  for  a  few  seconds,  run  for  52  seconds  at  a  time 
step  of  0.01  seconds,  and  sampled  every  .05  seconds  till  1024  samples 

of  input  to  the  man  and  output  from  the  plant  were  stored.  A  fast 
fourier  transform  was  taken  of  the  input  and  output.  The  open  loop 

gain  and  phase  were  calculated  from  zero  to  ten  radians  per  second  at 

intervals  of  one  radian  per  second  by  adding  the  complex  numbers  in 

every  eight  cells  and  then  dividing  the  absolute  values  to  get  gain 

and  subtracting  the  angles  to  get  phase.  The  computer  program  plots 
the  resulting  gain  and  phase  and  it  calculates  the  cost.  It  selects 
new  values  for  the  human  transfer  function  and  then  repeats  the 

simulation  and  cost  calculations  until  a  minimum  is  established. 

The  human  transfer  function  optimization  algorithm  worked  fine 

with  a  white  noise  input.  There  was  enough  power  at  every  frequency 
of  interest  to  make  good  calculations  of  gain  and  phase.  The  gain 

calculation  becomes  noisy  when  the  input  (the  denominator)  gets  near 
zero.  The  phase  calculation  gets  noisy  when  either  the  input  or 

output  power  gets  too  small  to  make  an  accurate  phase  measurement. 
The  algorithm  requires  smooth  monotonically-decreasing  gain  and  phase 
for  at  least  one  frequency  band  beyond  180  degree  phase  lag.  Such 
data  were  obtained  with  a  white  noise  input  with  and  without 
nonlinearities  in  the  loop.  Unfortunately  real  targets  do  not  present 
a  white  noise  tracking  spectrum. 

A  realistic  target  motion  was  constructed  from  the  following 
considerations: 

a.  The  algorithm  wants  as  much  power  as  possible  and  so  does  the 
maneuvering  target-,  therefore  a  course  made  of  segments  of  0.2g  turns 
was  used. 

b.  The  target  wants  to  move  forward  rather  than  go  in  a  circle; 
therefore  the  turns  were  limited  to  plus  and  minus  45  degrees  from  the 
line  of  sight  between  the  target  and  tracker. 
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c.  The  algorithm  wants  power  in  each  one-radian-per-second 
frequency  band,  therefore  the  turning  radius  and  speed  were  selected 
to  produce  a  fundamental  frequency  at  0.5  radian  per  second  so  that 

the  harmonics  would  fall  at  1.5,  2.5,  3.5 . radians  per  second. 

The  relationships  for  radial  acceleration  were  used. 

radial  acceleration  =  rw^  =>  v^/r  =  0.2g 
r  =  radius  a  8  meters 

w  a  angular  rate  *  0.5  radians  per  second 

g  *  acceleration  of  gravity  »  10  meters  per  second  squared 

d.  The  range  to  target  was  set  at  four  kilometers  to  reduce  the 
angular  tracking  rates  co  the  low  rates  of  interest. 

Maximum  rate  »  4  m/s  ♦  (sin  45)/4  km  *  0.7  mrad/s 

These  considerations  resulted  in  a  course  which  was  roughly  sinusoidal 
at  0.5  radians  per  second.  The  abrupt  changes  in  radial  acceleration 
every  90  degrees  of  turn  gave  strong  enough  harmonics  to  allow  the 
algorithm  to  calculate  gain  and  phase  when  there  were  no  nonlinear¬ 
ities  in  the  loop,  although  it  did  require  double  precision  in  the 
calculations.  When  nonlinearities  were  introduced  the  gain  and  phase 
curves  became  noisy  and  the  optimization  algorithm  would  not  work. 

Additional  power  was  required  at  both  the  input  and  output  at 
frequencies  of  interest  (1.5  through  8.5)  to  make  the  optimization 
algorithm  work  properly.  Fortunately  the  addition  of  noise  power  is 
justified  as  the  so-called  remnant  term  of  the  human  transfer  function 
(Yq) ,  Although  it  can  be  added  either  before  or  after  the  linear 
portion  of  Yq,  here  it  is  added  before  Yq  to  enhance  the 
optimization  algorithm  operation,  but  after  the  rms  calculation  to 
avoid  improperly  affecting  it.  The  appropriate  amount  of  noise  was 
calculated  by  the  following  steps: 

a.  Sheridan  page  241  shows  the  noise  power  to  be  20  percent  of 
the  total  power  at  the  output  of  the  man.  Page  242  shows  it  to  be 
uniform  with  frequency. 

b.  The  Yq  can  be  approximated  by  a  pure  gain  for  power 
calculations,  because  the  transportation  delay  does  not  affect  power 
and  the  lead- lag  terms  are  very  small. 

c.  Sample  trials  have  shown  the  rms  error  to  be  approximately 
0.3  milliradians . 

d.  There  are  approximately  eight  bands  of  interest. 


e.  Therefore  power  from  a  sinusoid  with  an  rms  amplitude  of  0.05 
milliradians  should  be  introduced  at  each  frequency  band. 

0.3  yo.2/8  -  0.05 


This  additional  power  helped  but  sooner  or  later  as  the  magnitude  of 
the  nonlinearities  was  increased  the  algorithm  would  become  too  noisy. 
There  are  still  a  couple  of  tricks  to  try,  i.e.,  longer  running  time  and 
extrapolation  of  the  phase  curve  to  180  degrees  lag  rather  than  inter¬ 
polation  as  was  done  here.  These  will  be  tried  in  the  Phase  II  effort.  The 
current  effort  was  finished  by  using  the  simple  expedient  of  minimizing  rms 
error  and  forgetting  about  the  phase  and  gain  calculations.  This  procedure 
raised  the  gain  until  the  system  went  unstable.  It  is  a  consistent  method 
but  it  is  probably  not  typical  of  human  operation. 


IV.  EFFECTS  OF  NONLINEARITIES  ON  THE  TURRET  RESPONSE 

Three  nonlinearities  were  added  one  at  time.  Deadspace  was  added  at  G3 
on  Figure  1.  It  corresponds  to  the  deadspace  in  a  gunner's  control  for  the 
first  couple  of  degrees  of  rotation.  Coulomb  friction  was  added  at  G4.  It 
corresponds  to  the  friction  on  the  turret  itself.  Backlash  was  applied  to 
the  output  at  H  and  it  can  also  be  applied  to  the  input  at  G3. 

Figure  2  through  Figure  6  show  the  effects  of  these  nonlinearities  on 
the  gain  and  phase  characteristics  of  the  turret,  i.e.,  from  G3  to  H. 
Figure  2  shows  the  turret  with  no  nonlinearities  for  a  comparison.  The 
turret  parameters  were  B  *  10,  KB  *  1.0,  and  J  *  1.0.  These  curves  were 
generated  by  using  a  single  sinusoid  by  itself  at  each  frequency.  The 
family  of  curves  in  Figure  3  through  Figure  6  represents  successive  doubling 


Figure  2.  Turret  fceeponse  Without  Honlineeritie*. 
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of  the  ratio  of  the  nonlinearity  to  the  input.  It  was  done  by  changing  the 
input  rather  than  the  turret.  This  procedure  is  equivalent  to  the  way  the 
test  data  would  be  collected  on  a  turret.  The  input  amplitudes  were  0.125, 
0.25,  0.5,  1,  2,  4  and  8  milliradians . 

The  phase  lag  shown  in  Figure  2  points  out  a  limitation  of  this 
methodology.  A  quick  calculation  would  predict  a  lag  of  135  degrees  at  10 
radians  per  second.  The  figure  shows  a  lag  of  150  degrees.  The  difference 
of  15  degrees  must  be  due  to  the  analysis  technique  which  uses  the  Zeta 
transform  and  the  Fast  Fourier  Transform.  The  time  step  used  with  the 
simulation  can  account  for  6  degrees  of  error  (0.01  seconds  x  10  radians  per 
second  x  60  degrees  per  radian).  The  rest  is  either  due  to  the  FFT  or  it  is 
unknown. 
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Figure  3.  Turret  Response  with  Coulomb  Friction. 


Figure  3  shows  the  effect  of  coulomb  friction  applied  to  G4.  The 
magnitude  was  0.1  pound-foot  applied  to  a  turret  of  one  slug-foot-squared 
polar  moment.  Hie  ratio  is  about  right  since  a  tank  has  a  22,000 
slug-feet-squared  polar  moment  and  about  2000  pound-feet  of  coulomb  friction 
referred  to  the  turret.  The  gain  decreased  as  the  input  was  decreased  but 
the  phase  lag  decreased  as  well. 
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FREQUENCY  IN  RRDIRN5  PER  SECOND 

Figure  4.  Turret  Response  With  Dead  Space  at  the  Control  Handle. 


Figure  4  shows  the  effect  of  deadspace  applied  at  G3.  The 
magnitude  of  the  deadspace  was  0.1  milliradians .  A  typical  turret 
might  have  0.04  radians  deadsDace  at  the  turret  control  handle.  The 
turret  gain  during  these  turret  response  runs  was  0.1  compared  to  a 
typical  turret  gain  of  0.02  radians  per  second  per  radian.  Obviously 
the  problem  will  require  new  coefficients  for  a  quantitative  analysis 
but  these  figures  show  the  trends. 
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FREQUENCY  IN  RR01RN5  PER  5EC0ND 


Figure  5.  Turret  Response  With  Backlash  at  the  Control  Handle. 


Figure  5  shows  the  effect  of  backlash  applied  at  the  turret 
control.  The  magnitude  of  the  backlash  was  0.04  milliradians . _  Once 
again  the  level  chosen  was  not  necessarily  representative  ot  real 
turrets,  however  it  does  show  the  relative  effects  of  backlash  on  gain 
and  phase.  Backlash  at  the  control  handle  will  cause  a  nhase  <g 
without  changing  the  gain  substantially. 


15 


^  «C  ,^*—V  1  .MM  .  *  ~  ,  - 


Figure  6.  Turret  Response  With  Backlash  at  the  Turret  Output. 


Figure  6  shows  the  effect  of  0.001  milliradian  of  backlash  on  the 
output.  Here  the  backlash  can  be  seen  to  have  a  greater  effect  at 
higher  frequencies  as  compared  to  backlash  on  the  input.  The 
reduction  in  amplitude  with  increased  frequency  at  the  output  causes 
this  effect. 

The  objective  of  presenting  these  figures  is  to  indicate  that  it 
will  be  possible  to  shape  the  gain-phase  characteristics  of  the  turret 
model.  This  will  be  done  when  the  test  data  from  the  tank  turret 
become  available. 
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V.  EFFECTS  OF  NONLINEARITIES  ON  LOOP  RESPONSE 


The  intent  at  this  point  was  to  calculate  the  rms  error  for  the 
closed-loop  system  when  tracking  the  target  course  developed  earlier.  The 
gain-phase  plots  for  the  open-loop  response  of  the  man  and  turret  were  also 
of  interest,  but  as  explained  earlier  the  gain-phase  plots  were  usable  only 
for  the  condition  with  low  levels  of  nonlinearities.  Figure  7  shows  these 
plots  for  a  condition  with  no  nonlinearities.  This  condition  had  a  phase 
margin  of  41  degrees,  a  gain  margin  of  4.8  db  and  an  rms  error  of  0.31 
milliradians .  The  cross-over  frequency  (0  db  gain)  was  3.0  radians  per 
second . 


Frequency*  radians  per  second 


Figure  7.  Gain-Phase  Plot  Without  Nonlinearities. 


When  the  nonlinearities  were  added  the  loop  was  optimized  for  minimum  • 

rms  error.  The  rms  error  for  no  nonlinearities  dropped  to  0.26  milliradians  ■ 

but  the  loop  was  not  nearly  as  stable.  The  phase  margin  was  only  13  degrees  j 

and  the  gain  margin  was  1.4  db.  The  growth  in  rms  error  with  increased  j 

levels  of  nonlinearities  is  shown  in  Figures  8,  9,  and  10.  The  deadspace  in  j 

Figure  8  is  at  the  control  handle.  The  turret  gain  was  changed  to  0.01  (BK  j 

=  0.10,  B  =  100)  for  these  runs.  The  coulomb  friction  in  Figure  9  was  j 

applied  to  G4.  The  backlash  in  Figure  10  was  applied  to  the  control  $» 

handle .  4 
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DEAD  SPACE,  RADIANS 


re  8.  Tracking  Error  with  Dead  Space  at  the  Control  Handle. 


9.  Tracking  Error  With  Coulomb  Friction  at  the  Turret  Output 


mtur*%Pftng-f*r  ^ 


Figure  10.  Tracking  Error  with  Backlash  at  the  Control  Handle. 


The  only  conclusion  that  can  be  reached  at  this  time  regarding  Figures 
8,  9,  and  10  is  that  the  level  of  nonlinearities  that  were  used  did  have  an 
influence  on  tracking  error.  It  remains  to  be  seen  if  these  are  the 
appropriate  levels.  The  turret  measurement  tests  will  determine  the 
appropriate  levels  to  use.  It  also  remains  to  be  proven  that  the  man-model 
used  here  is  appropriate  for  this  task.  The  human  tracking  tests  will 
determine  that. 


VI.  SUMMARY 

A  control  loop  with  a  man-model  and  with  provision  for  nonlinearities 
was  developed.  An  optimization  algorithm  for  the  adaptive  man-model 
worked  well  for  low  levels  of  nonlinearities,  but  it  had  to  be  simplified  to 
work  for  high  levels  of  nonlinearities.  Nonlinearities  oere  shown  to 

influence  tracking  error. 

The  next  phase  of  this  effort  will  have  the  benefit  of  quantitative 
descriptions  of  the  turret  response.  The  turret  will  be  simulated  in  more 
detail  and  the  correct  parameter  values  will  be  used  for  nominal  conditions. 
Another  attempt  will  be  made  to  improve  the  adaptive  man-model  to  work  wi^h 
the  appropriate  nonlinearities. 
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PROGRAM  CUNTRL 

THIS  PROGRAM  SIMULATES  A  TANK  GUNNER  TRACKING  A  TARGET. 
ON  DEMAND  IT  WILL  ADJUST  THE  PARAMETERS  OF  THE  MAN 
TO  GIVE  THE  LOWEST  COST  OR  ERROR. 

IT  WILL  ALSO  PRODUCE  BODE  PLOTS  ON  PRINTER  AND  IN  A 
FILE  FOR  USE  BY  A  PLOTTER. 

THE  PROGRAM  IS  INTERACTIVE,  AND  PROMPTS  ALL  INPUTS. 

A  CARRIAGE  RETURN  IS  SUFFICIENT  FOR  AN  ANSWER  OF 
ZERO  OR  NO. 

TWO  PROMPTS  REQUIRE  MULTIPLE  INPUTS 
ON  ONE  LINE,  SEPERATED  BY  COMMAS. 

FIRST  SET: 

TRAILING  ZERO  VALUES  MAY  BE  IGNORED, 

COMMAS  ARE  SUFFICIENT  FOR  NONTRAILING  ZERO  FIELDS. 
TAU  TRANSPORT  DELAY  IN  SECONDS 

B  VISCOUS  FRICTION 

BK  PLANT  GAIN 

BKLSH  BACKLASH  AT  OUTPUT  IN  RADIANS 

CF  COULOMB  FRICTION 

DEDSPC  DEADSPACE  AT  INPUT  IN  RADIANS 


SECOND  SET: 
TI 
K 
TL 


INTEGRATION  TIME  IN  SECONDS 

GAIN  OF  NAN 

LEAD  TIME  IN  SECONDS 


INSTEAD  OF  THE  SECOND  SET,  DEFAULTS  OF 
TI  =  .01,  K  =  2.5  *  B,  AND  TI  =  TI  +  1/B 
MAY  BE  CALLED  BY  A  CARRIAGE  RETURN. 

THE  PLANT  IS  NORMALIZED  TO  A  MASS  OF  1. 

SUBROUTINES  IN  T1IE  PACKAGE: 

THE  MAIN  PART  OF  THE  PROGRAM  DOES  AI.L  THE  INTERACTIVE 
CONVERSATION  AND  CALLS  PLANT©,  AUTO,  FUN,  AND  MACHINE. 

SUBROUTINE  AUTO 

SETS  UP  THE  AUTOMATIC  OPTIMIZATION 
CALLS  FUN  AND  FNMIN. 

THE  PARAMETER  ACRCY  IN  AUTO  TELLS  FNMIN 
THE  PRECISION  DESIRED. 

THE  AUTOMATIC  MINIMIZATION  ALSO  TERMINATES  IF  TI 
BECOMES  LESS  THAN  .005 


SUBROUTINE  FNMIN 

DOES  THE  AUTOMATIC  OPTIMIZATION 

CALLS  FUN 

IF  AUTOMATIC  MINIMIZATION  IS  CHOSEN 
FNMIN  SYSTEMAT I CALY  VARIES 

X(  1  >  =  1  /  TI ,  X(  2  >  =  K,  AND  X(  3  >  =  TL 
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o  TO  MINIMIZE  THE  COST  RETURNED  BY  FUNCTION  FUN. 

c  X(  1  )  IS  INVERTED  TO  AVOID  NEGATIVE  VALUES, 

c 

c  FUNCTION  FUN 

o  MODELS  THE  CONTROL  LOOP 

o  COMPUTES  THE  COST  OF  TRACKING 

o  CALLS  MAN©,  PLANT 1 ,  STATS0 ,  MANI ,  TGTS,  STATS 1 ,  PLANT2, 

c  STATS2 ,  STATSW,  AND  FFT. 

o  THE  PARAMETER  NN  IS  USED  AS  A  FLAG  IN  FUN: 

o  NN  >  0,  AUTOMATIC  REDUCTION,  NO  BODE  PLOT, 

o  NN  3  0,  BODE  PLOT  OF  MAN-MACHINE  SYSTEM 

o  USING  THE  MODEL  TARGET  AS  INPUT, 

o  NN  3  -1,  BODE  PLOT  OF  MACHINE  WITH  MODEL  TARGET 

o  AS  INPUT  TO  MAN, 

o  NOTE  THAT  IN  THIS  CASE  THE  TARGET  IS  FILTERED  THRU 

o  THE  MAN  AND  THE  PLOT  IS  THEREFORE  AN  IMPLICIT  FUNCTION 

o  OF  THE  MAN. 

o  NN  <  -1,  BODE  PLOT  OF  MACHINE  WITH  SINE  WAVE  INPUT. 

c 

o  SUBROUTINE  MACHINE 

o  MAKES  BODE  PLOTS  OF  THE  MACHINE 

o  CALLS  FUN  AND  PLOT, 

c 

c  SUBROUTINE  PLOT 

c  PRODUCES  A  PRETTY  BODE  PLOT  IN  A  FILE  READY  FOR  PLOTTING 

o  PLOT  ASSUMES  THE  PLOTTING  PACKAGE 

c  TIG  (  TERMINAL  INDEPENDENT  GRAPHICS  )  WHICH  WAS 

c  Will  TEN  IN  C  AND  REQUIRES  A  C  COMPILER. 

c 

c  SUBROUTINE  FFT 

<-  FAST  FOURIER  TRANSFORM 

o 

c  THE  MAN : 

c 

c  SUBROUTINE  NANO 

c  INITIALIZES  THE  MAN 

c  CALLS  TGTAO 

c 

o  SUBROUTINE  MANI 

c  THE  HAN’S  PART  OF  THE  CONTROL  LOOP 

c  CALLS  TCTA1 

c 

c  THE  PLANT : 

o 

o  SUBROUTINE  PLANTO 

c  INITIALIZES  THE  PLANT 

c 

c  SUBROUTINE  FLANT1 

c  HZSETS  PLANT  AT  START  OF  EACH  RUN 

c 

c  SUBROUTINE  PLANT2 

c  THE  MACHINERY’ S  PART  OF  THE  CONTROL  LOOOP 

o 

c  THE  STATS: 
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o 

o  SUBROUTINE  STATS 0 

c  INITIALIZES  STATISTICS  PACKAGE 

o 

c  SUBROUTINE  STATS 1 

c  COLLECTS  THE  DATA 

c 

o  SUBROUTINE  STATS2 

o  FINDS  RMS  ERROR,  MEAN  INPUT  OFFSET,  ETC 

o 

o  SUBROUTINE  STATSW 

o  WRITES  THE  STATS 

o 

o  THE  TARGET: 

o 

c  SUBROUTINE  TGTA0 

o  SETS  UP  TARGET 

o 

c  SUBROUTINE  TGTA1 

o  RUNS  TARGET 

o 

o  SUBROUTINE  TGTS 

o  SINE  WAVES  FOR  BODE  PLOTS 

o 

o  THE  PROGRAM  WAS  WRITEN  TO  BE  RUN  ON  A  PDP  11/70 

o  USING  THE  CULC  F4P  COMPILER  AND  THE  UNIX  OPERAT’iNG 

o  SYSTEM. 

o  THE  PROCRAM  IS  IMPLICIT  DOUBLE  PRECISION, 

c  AND  USES  COMPLEX  ARITHMATIC. 

c  GENERIC  NAMES  HAVE  BEEN  USED  FOR  FUNCTIONS,  IE  ABS, 

c  AND  THE  F4P  COMPILER  SELECTS  THE  PROPER  FUNCTION, 

c  IE  ABS.  CABS,  DABS,  ETC. 

c 
c 

common  b,  delta,  del  tot, 

+  Jmax,  Jmaxt,  Jmod,  k,  nc ,  nt ,  pi,  tmnx,  uni  In 
common  /  x  /  x(  10  > 
external  fun 
real  *  0  JJ,  kk 

write  (  6 ,  5  ) 

5  format  (  *  NOW  TYPE  VALUES  FOR  TAU,  B,  BK,  BACKLASH,  ’  / 

+  ’  COULOMB  FRICTION,  AND  DEAD  SPACE  ’  ) 

accept  10,  tau,  b,  bk,  bklsh,  cf,  dedspc 
10  format  (  7fl0.0  ) 

uni  In  =  abs(  bklsh  )  +  abs(  cf  )  +  abs(  dedspc  ) 

k  =  10 

Jmod  =  5 

Jmax  =  2  **  k 

jmaxt  =  jmax  *  jmod 

tetep  =  2. 

pi  =  3.1415926535 

nc  =  8 

mi  =  0 
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fno  =  no 

tmax  s  fno  #2.  *  pi 

delta  »  tmax  /  float!  Joax  ) 

deltat  ‘  tmax  /  float!  Jmaxt  ) 

write  !  6,  15  )  tnu,  delta,  b,  bk ,  bklsh,  of,  dedapo,  Jmax 
15  format  (  /  ’  TAU  DELTA  D  BK  BKL8H  CLFR  DEDSP  JMAX  *  / 

+  2f6.3,  f6.1,  f 6 .3,  3pf6.2,  0pf6.3,  3pf6.2,  15  /  ) 

nt  =  tau  /  doltat  +  .5 

call  plantO!  b,  bk,  bklsh ,  bklshl,  of,  dedapo,  deltat  ) 
write  !  6,  20  ) 

20  format  !  ’  DO  YOU  WANT  TO  LOOK  AT  THE  MACHINERY?  *  / 

+  ’  TO  ANSWER,  TYPE  1  FOR  YES,  0  FOR  NO.  '  ) 

aooept  30,  m 
write  !  6,  25  ) 

25  format  !  /  •  DO  YOU  WANT  AUTOMATIC  REDUCTION  ?  *  / 

+  *  TYPE  1  FOR  YES,  OR  0  FOR  NO  ’  ) 

aooept  30,  n 
30  format  !  110  ) 

35  write  !  6,  40  > 

40  format  !  /  '  NOW  TYPE  VALUES  FOR  TI ,  K,  AND  TL  ’  ) 

accept  10,  tl,  kk,  tl 
If  !  tl  .gt.  0.  )  go  to  45 
tl  =  .01 
kk  =  2.5  *  b 
tl  =  tl  +  1 .  /  b 
45  x!  1  )  a  1  .  /  tl 

x(  2  )  =  kk 

x(  3  )  =  tl 

If  l  n  .gt .  0  )  call  auto 

cost  =  fun!  x,  nn  ) 

If  (  m  .eq.  1  )  call  machine 

If  (  n  .le.  0  )  go  to  35 

stop 

end 

subroutine  auto 

Implicit  doublo  precision  (  a-h,  o-z  ) 
common  b,  delta,  deltat, 

+  Jmax,  jmaxt,  Jmod,  kj,  nc ,  nt,  pi,  tmax,  uni  In 
dimension  ops!  10  ) 
common  /  x  /  x!  10  ) 
real  *  0  kk 
external  fun 

k  =  1 
nn  =  0 

cost  =  fun!  x,  nn  ) 
tl  =  1 .  /  x<  1  ) 
kk  =  x  (  2  ) 
tl  =  x(  3  ) 

write  (  6,  15  )  cost,  tl,  kk,  tl 
15  format  (  /  ’  THE  INITIAL  COST  WAS  ’ ,  3pf8.3  / 
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+  ’  FOR  TI  =  * ,  0pf7.4,  ’  KK  =  ’ ,  f5.1 ,  ’ ,  TL  =  ’ ,  f7.4  // 

+  ’  TI  KK  TL  RMS  PHASMR  GAINMR  COST  SLOPE’  /  ) 

acrcy  =  .05 
e  =  5 . 
n  =  3 

do  5  1  -  1 ,  n 

eps!  i  )  =  acrcy  *  x!  1  ) 
call  fnmin!  n,  x ,  cost,  fan,  e ,  tps ,  k  ) 
tl  =  1 .  /  x!  1  ) 
kk  =  x<  2  ) 
tl  =  x(  3  ) 

write  (  6,  10  )  cost,  tl,  kk,  tl 

format  !  /  ’  THE  MINIMUM  COST  WAS  ’  3pf6.3  / 

+  ’  FOR  TI  3  ’  ,  0pf7 .4,  * ,  KK  3  ' ,  f5. 1 ,  ’ ,  TL  =  ’ ,  f7.4  ) 

return 

end 

function  fun(  x,  nn  ) 

Implicit  double  precision  (  a-h,  o-z  ) 
common  b,  delta,  deltat, 

+  Jmax,  Jmaxt ,  Jmod,  k,  no,  nt ,  pi,  tmax,  uni  In 
common  /  gp  /  gnan!  10,  7  ),  phaas(  10,  7  ) 
dimension  gain!  81  ),  phase!  81  ) 
dimension  gl  !  100  ) 

dimension  frlnC  1024  ),  flin(  1024  ),  frout(  1024  ), 

+  flout!  1024  ) 
dimension  x!  10  ) 

complex  oln,  cout,  orln,  crout,  sorln,  sorout,  oxx,  cyy 
real#4  filnl,  flout!,  frlnl ,  froutl 
data  cos tmn  /  1  .  / 


If 

! 

X  ( 

1 

)  .le. 

0. 

.or.  x!  2  ) 

.  le . 

o. 

•or.  x!  3  ) 

.  le. 

o. 

) 

cos  t 

= 

1 

.  e+2 

If 

! 

x( 

1 

)  .le. 

0. 

•or.  x!  2  ) 

.  le « 

o. 

. or .  x!  3  ) 

.  le  . 

0. 

) 

+  go  to  900 

call  manO!  deltat,  nt ,  x  ) 
call  plan' M  h  ) 
points  =  jmax 

If  !  nn  .oq.  0  )  write  (  6,  55  ) 

format  !  /  27x,  4h-1.0,  6x,  4h-0.5,  7x,  3h0.0,  7x,  3h0.5, 

+  ’  1.0’  / 

+  ’  SEC  REF  IN  OUT  +’ ,  4!  9x,  lh+  )  ) 

If  !  nn  .cq.  -1  )  write  !  6,  56  ) 

format  !  27x,  ’-250  -125.  0.0  125.’ 

+  ’  250’  / 

+  ’  SEC  REF  IN  OUT  +’  ,  4!  9x,  il»+  )  ) 

call  statsO!  delta,  Jmax  ) 

Jem  =  2 
Jmx  =  Jmaxt 
Jcl  =  Jem 

If  !  nn  .go.  -1  )  Jcl  =  1 

NOW  MODEL  THE  CONTROL  LOOP 
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do  101  Jot  =  Jol,  Jem 
Jm  =  0 

if  (  jet  .  eq.  Jem  )  Jmx  =  jmaxt 
do  100  j  3  1,  Jmx 

If  (  nn  . go .  -  1  )  call  raanl (  b,  J,  nt ,  dlffdt,  dlfft,  g3  ) 

If  (  nn  .It.  -1  )  call  tgts!  J,  g3  ) 
if  (  jet  .It.  Jem  )  go  to  90 

If  (  mod(  J,  Jmod  )  ,ne.  0  )  go  to  90 
Jm  3  Jm  +  1 

if  (  nn  .ge.  0  )  frln(  Jm  )  3  dlffdt 

If  (  nn  .It.  0  )  frln(  Jm  )  3  g3 

frout!  Jm  )  3  h 
flln(  Jm  )  3  0. 
flout(  Jm  )  3  0. 

If  (  nn  .go.  -1  ) 

+  call  statsl!  dlfft,  frin(  Jm  ),  frout(  Jm  ),  Jm  ) 

90  call  plant2(  g3,  b  ) 

100  continue 

101  continue 

If  (  nn  .go.  -1  ) 

+  call  stats2!  avcln,  aveout,  rms  ) 
do  120  J  3  1 ,  Jrnax 
tj  =  J 

frlnC  J  )  3  frln!  J  )  -  aveln 
frout(  j  )  3  frout!  J  )  -  aveout 

120  continuo 

121  continue 

If  (  nn  .gt.  0  .or.  nn  .It.  -1  )  go  to  170 
do  150  J  =  1,  199,  3 
t  3  float!  J  )  *  delta 

If  (  nn  . eq.  0  )  n  3  frln(  J  )  *  20000. 

if  (  nn  .It.  0  )  n=  frin(  J  )  *80. 

ref  3  frln(  J  )  +  frout(  j  ) 

If  (  n  .gt.  -21  ) 

+  write  (  6,  140  )  t,  ref,  frln(  j  ),  frout(  J  ) 

If  (  n  .le.  -21  )  write  <  6,  142  )  t,  ref,  frin!  J  ),  front!  J  ) 
150  continue 

140  format  !  lx,  f6.2,  3p3f7.2,  <n+21>x,  lb*  > 

141  format  (  lx,  4el0.3  ) 

142  format  !  f7.2,  3p3f7.2  ) 

170  if  !  nn  .eq.  0  )  coll  statsw 
175  call  fft  !  frln,  flin,  k  ) 

call  fft  !  frout,  flout,  k  ) 
nav  3  nc 

if  !  nn  .It.  -1  )  nav  =  1 

Irads  3  1 
1 rndm  3  42 

If  !  nn  .It.  -1  )  Iradm  3  1  -  nn 
If  !  nn  .It.  -1  )  Irnds  3  iradm 
do  500  Irad  3  Irads,  Iradm 
serin  3  0. 
soront  3  0. 
do  400  J  3  1 ,  nav 

1  3  nav  *  !  1  rad  -  1  )  +  j 
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frlnl  =  frln(  1  ) 

f i Ini  3  f 1 ln<  1  ) 

cln  =  oraplx!  frlnl,  fllnl  ) 

froutl  =  frout(  1  ) 

fioutl  =  flout!  1  ) 

cout  =  cmplxt  froutl,  flontl  ) 

crln  =  cln  *  conjg!  cln  ) 

crout  =  oout  *  oonjg!  cln  ) 

If  <  abs  l  orln  )  .le.  1 .e-10 
+  .or.  nn  .It.  -1  .or.  nn  .gt.  0 
+  .or.  1  .gt.  80  )  go  to  391 
oxx  s  orout  /  orln 
cyy  =  nbs!  oxx  ) 

gan  3  20.  *  logl0(  renl(  cyy  )  ) 

faz  3  atnn2(  alinag!  orout  ),  real!  crout  )  ) 

faz  =  180.  *  faz  /  pi 

enl  3  real!  orln  ) 

onl  3  sqrt!  onl  ) 

oyy  3  abs!  oout  ) 

eno  3  real!  oyy  ) 

If  !  nbs(  ono  )  .gt.  1 .©-10  )  write  t  6,  390  ) 

+  1,  frln!  1  ),  flint  1  ),  frout!  1  ),  flout!  1  ), 

+  gnu,  fnz,  onl,  ono 

390  formnt  !  IS,  4f0.4,  2f8.2,  2f8.4  ) 

391  contl nuo 
ho rln  3  serin  +  crln 
scrout  3  scrout  +  crout 

cont 1 nuo 

If  t  scrout  .  ©q .  0.  )  write  t  6,  450  ) 

If  t  obs!  scrout  )  .gt.  1 .e+1©  )wrlte  t  6,  431  ) 

1 .  /  x<  1  ) ,  x(  2  ) ,  x(  3  ) 
format  t  /  f7.4,  f7.2,  f7.4,  *  UNSTABLE  *  /  ) 
fonnut  t  /  •  INPUT  TOO  SMALL  ’  ✓  ) 

If  t  scrout  .oq.  0.  .or.  nbn!  scrout  )  .gt.  l.e+10  ) 
cost  3  1 . o+l 

If  !  scrout  .oq.  0.  .or.  nbst  scrout  )  .gt,  l.o+lO  ) 
go  to  900 
exx  3  0. 

If  t  nbs!  scr\n  )  .gt .  l,o-30  )  oxx  3  scrout  /  serin 
cyy  3  nbs!  exx  ) 
gain!  1  rad  )  3  real!  cyy  ) 

phase!  Irnd  )  3  atan2!  nlmagt  scrout  ),  real!  scrout  >  ) 
If  t  phase!  Irnd  >  . gt.  0.  ) 
phase!  Irnd  )  3  phase!  Irnd  )  -  2.  *  pi 
cont 1 nuc 

If  t  nn  .gt.  0  .or.  nn  .It.  -nc  )  go  to  540 

If  !  nn  .eq.  0  )  Idbmn  3  -10 

If  !  nn  .eq.  -I  )  Idbmu  3  -60 

If  t  nn  .It.  -1  )  Idbmn  3  -70 

Idbmx  3  Idbmn  +  30 

write  t  6,  510  )  t  Idb,  Idb  3  Idbmn,  Idbmx,  10  ) 
format  !  /  20x,  ’PHASE  SHIFT’,  20x,  71iDB  GAIN  / 

’  FREQ  PH  GAIN  -180’,  15x,  ’-90*. 

14,  3!  7x,  13  )  ) 


400 

+ 

431 

450 

+ 

+ 


+ 

500 


510 

+ 

+ 
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writ©  (  6,  530  ) 

530  format  (  lx,  lh+,  4<  8x,  lh+  ),  3x,  lh+,  3(  9x,  Ih+  )  ) 

540  xx  =  180.  /  pi 

radf  =  2.  *  pi  #  float!  nnv  )  /  tmax 
npm  =  52 

If  (  nn  .It.  0  )  npm  =  n^ra  -  Idbmn  -  10 
11  a  1 

lm  =  40 

if  (  nn  .It.  -1  )  11  «  1  -nn 
If  (  nn  .It.  -1  )  lm  =  1  -  nn 
do  600  1  s  11,  lm 

If  <  gain!  1  )  .le.  1.  .or.  galn(  1+1  )  .gt.  1.  ) 

+  go  to  590 

phasmr  =  180.  +  xx  *  (  phase!  1  )  - 
+  (  phase!  1  )  -  phase!  1+1  )  ) 

+  *(galn(l)-l.  ) 

+  /  (gain!  1  )-galn(  1+1  )  )  ) 

If  (  l  .eq.  1  ) 

+  slope  »  20.  *  loglO!  gain!  1  )  /  gain!  l+l  )  ) 

+  /  loglO!  1.  /  3.  ) 

If  (  1  .gt.  1  ) 

+  slope  =  20.  *  loglO!  gain!  l-l  )  /  gain!  1+2  )  ) 

+  /  loglO!  (  float!  1-1  )  -  .5  ) 

+  /  (  float!  1+1  )  +  .3  )  ) 

590  n  =  phase!  1  )  *  xx 

up  =  n  /  5 

gnn  =  20.  *  loglO!  gain!  1  >  ) 
ra  =  gnn 

rf  =  radf  *  (  float!  1)  -  .5  ) 

If  (  nn  .It.  -1  )  rf  =  rf  -  .5  /’  float!  no  ) 
nr f  -  rf  +  .001 

If  (  nn  .oq.  -no  )  111  *111+1 

If  (  nn  .go.  -1  .or.  nrf  .It.  1  .or.  Ill  .It.  1  )  go  to  595 
If  (  nn  .It.  -1  )  gaan!  nrf,  ill  )  =  gnn 
If  (  nn  .It.  -1  )  phnns!  nrf,  111  )  =  phase!  1  )  *  xx 
595  continue 

If  (  nn  .lo.  0  .and.  np  .gt.  -55  .and.  m  .gt.  2-npir  ) 

+  write  (  6,  610  )  rf ,  n,  gnn 

If  (  nn  .lo.  0  .and.  (  np  .io.  -55  .or.  m  .lo.  2-npm  ) 

+  .and.  gnn  .gt.  -100.  ) 

+  write  (  6,  611  )  rf,  n,  gan 

If  (  phase!  I  )  .gt.  -pi  .and.  phase!  1+1  )  .It.  -pi  ) 

+  gultunr  =  -  20.  *  loglO!  gain!  1  ) 

+  -  (  gain!  I  )  -  gain!  l+l  )  ) 

+  *  (  plinso!  1  )  +  pi  ) 

+  /  (  phase!  1  )  -  phase!  l+l  )  )  ) 

If  (  nn  .It.  - 1  )  go  to  600 

If  (  phase!  1  )  .It.  -pi  .and.  gain!  1  )  .It.  1.  )  go  to  700 
600  continue 

610  format  (  lx,  f5.2,  15,  f5.l,  t<np+56>,  lh*,  t<m+npm>,  lh+  ) 

611  format  (  lx,  f5.2,  15,  fS.l  ) 

620  f o rmn t  (  lx,  flO.l,  110,  fl0.3  ) 

If  (  nn  . go .  - 1  )  write  (  6 ,  630  ) 

630  format  (  ’  FELL  THRU  600  LOOP  ’  ) 
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700  coot  =  rmo 
go  to  701 

If  (  phasrar  . ne.  0.  ) 

+  cost  =  rmo  /  sln(  phasrar  *  pi  /  90.  ) 

If  (  gnlnmr  .It.  6.  .and.  gainrar  .gt.  0.  ) 

+  cost  =  coot  /  sln(  galnmr  /  6.  *  pi  /  2.  ) 

701  continue 

If  (  nn  .  eq.  0  ) 

+  write  (  6,  750  )  phasmr,  galrunr,  oost,  slope 
750  format  (  /  *  THE  PHASE  MARGIN  IS  ’  ,  f5.0,  *  DEGREES  ’  / 

+  *  THE  GAIN  MARGIN  IS  f 6. 1 ,  *  DB.  ’  / 

+  *  THE  COST  IS  * ,  3pf 12 .3  / 

+  *  THE  SLOPE  IS  ’  ,  0pf6 . 1 ,  •  DB  PER  DECADE.  •  /  ) 

tl  »  1.  ✓  x(  1  ) 

if  (  nn  .gt.  0  .and.  coat  .It.  1.  )  write  (  6,  800  ) 

+  1 1 ,  x(  2  ) ,  x<  3  )  ,  rmo,  phasmr,  galrunr,  oost,  slope 

800  format  (  f7.4,  f7.2,  f7.4,  3pf7.3,  0p2f7.3,  3pf  3,  0pf7.2  ) 
If  (  tl  .It.  .005  .and.  cost  .It.  costmn  )  cost  5  -  cost 
oostmn  =  mln(  cost,  costmn  ) 
fun  -  oost 
900  return 


and 

subroutine  machine 

Implicit  double  proolslon  (  n-lr,  o-z  ) 
common  b,  delta,  del  tut, 

Jmnx ,  jmaxt,  Jmod ,  KJ ,  nc ,  nt ,  pi,  tmax,  nr.lln 
oommon  /  gp  /  gaunt  10,  7  ),  phanst  10,  7  ) 
common  /  mark  /  mnrk(  8  ) 
common  /  tgts  /  a,  w 
oommon  /  x  /  x<  10) 
byte  or 

data  mark  /  ’+’  ,  » *’  ,  ’ x’ ,  ’o*  ,  ’s’  ,  ' ,  *  .’  ,  ’  n‘  / 
or  =  "015 


wrlto  (  6,  5  ) 

5  format  (  /  *  NOW  LOOK  AT  MACHINE  ONLY  ’  /  ) 

n  =  -1 

cost  =  fun(  x,  n  ) 
write  (  6,  10  ) 

10  format  (  /  ’  BODE  PLOTS  FOR  SINUSOIDAL  DRIVING  FUNCTIONS  ’  /  ) 

do  200  k  =  1,  7 

a  =  . I  /  2  **  (  k-1  ) 
wrlto  (  6,  20  )  a 

20  format  (  /  lx,  3pf6.2,  ’  MILLIRADIANS  PER  SECOND  AMPLITUDE.  ’  ) 

km  =  k 

do  100  l  =  1 ,  10 

w=  del  tat  *  float(  1  ) 

n  =  -  nc  *  1 

oost  =  fun(  x,  n  ) 

if  (  cost  .gt.  1.  )  go  to  300 

100  continue 
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if  (  unlin  .eq.  0  )  go  to  300 

200  continue 

300  write  (  6,  320  ) 

320  format  (  //  ’  CONSOLIDATED  PLOT  ’  // 

+  20x,  ’  PHASE  SHIFT  ’ ,  20x,  ’  DB  GAIN  * 

+  /  *  FREQ  -225  -1B0  -135  -90  ’ 

+  ’  -70  -60  -30  -40*  / 

+  10x,  lh+,  3(  Bx,  lh+  ),  3x,  lh+,  3(9x,  lh+  )  ) 

do  400  1=1,10 
write  (  6,  330  )  1 
330  format  (  15,  •  ) 

itpm  =  uwrlte  (  1,  or,  1  ) 

do  350  k  =  l ,  km 
m  =  gaan(  1 ,  k  ) 
np  =  phaas(  1,  k  )  /  5. 
if  (  m  .gt.  -100  .and.  np  .ge,  -54 
+  . and .  m  . 1 t .  0  . and .  np  .It.  0  ) 


+ 

write  (  6,  340  )  mark! 

k  ),  mark!  k  ) 

340 

format  (  t<np+56>,  al , 

t <m+l 12> ,  al,  8 

Itpm  =  uwr 1 1 o (  1 ,  c r 

1  ) 

350 

cont Inue 

write  (  6,  360  ) 

360 

format  (  lx  ) 

400 

cont 1 nue 

call  plot(  km  ) 
return 


end 

subroutine  plot(  n  ) 

implicit  double  precision  (  a-h,  o-z  ) 
common  /  gp  /  gaan(  10,  7  ),  phaas(  10,  7  ) 
common  /  murk  /  mnrk(  8  ) 

dimension  lxp(  10  >,  lyp(  10  ),  lyg(  10  ) 
byte  dov(  4  ),  fllot  10  ) 

external  ffaxls,  fflino,  ffnoworlgln,  ff out put 
data  dev  /  ’v’ ,  ’ t’ ,  ’o’ ,  0  / 

data  file  /  ’d’,  ’a*.  *  t ’ ,  ’a’,  ’p’,  ’  l  ’  ,  ’o’,  ’t’, 
+  2*0/ 


coll  cnllci  ffoutput,  file,  0  ) 
coll  cnllcl  ffnoworigin,  1000,  1000  ) 
call  cal lo(  ffaxls,  ’  FREQUENCY  IN  RADIANS  PFJ\  SECOND 
+  0,  0,  6000,  0.0,  0,  0.,  1.,  600  ) 

call  cnllc(  ffaxls,  ’  PHASE  IN  RADIANS  ==  GAIN  IN  I)R 
+  0,  0,  -6000,  90.,  0,  -200.,  20.,  600  ) 

Jmx  =  10 
do  100  1  =  1  ,  n 
do  90  J  =  1  ,  10 

lyg<  J  )  =  gnan(  j,  1  )  *  30.  +  6000. 

lyp(  j  )  =  phnas(  J,  1  )  *  30.  +  6000. 

1 xp (  j  )  =  J  *  600 
cont 1 nue 


90 


Jan  2  10:41 


call  callc(  ffllne,  Ixp,  iyp,  10,  nark(  i  ),  1  ) 
call  callc(  ffllne,  Ixp,  lyg,  10,  roark(  1  ),  1  ) 
100  continue 
return 

end 


aubroutlne  tgtaO 

Implicit  double  precision  (  a-h,  o-z  ) 
common  b,  delta,  deltat, 

+  Jmax,  Jmaxt ,  jmod,  k,  no,  nt,  pi,  tmax,  uni  In 
ooramon  /  tgtn  /  dt,  h,  hdt,  hllm,  ood,  pos,  vdt 

dlst0  3  4000. 
dt  3  del  tut 
hllm  3  pi  /  B, 
pos  3  0. 
h  3  hllm 
ve 1  3  10. 
hdot  3  vol  /  30. 
hdt  3  hdot  *  deltnt 
ood  3  1 .  /  dls tO 
vdt  3  vol  *  deltnt 
return 
o 

end 

o 

subroutine  tgtnl (  J,  dither,  g  ) 
o 

1'ipllclt  double  precision  <  a-h,  o~z  ) 
common  /  tgtn  /  dt,  h,  hdt,  hllm,  ood,  pos,  vdt 
o 

If  (  nbsl  h  )  .go.  hllm  .and.  hfchdt  .gt.  0.  ) 

+  hdt  3  -  hdt 
h  3  h  +  hdt 

pos  =  pos  +  sln(  h  )  *  vdt 
dt  f  3  dt  *  float (  J  ) 
angle  3  ood  *  pos 

dither  3  ood  *  (  .1  *  sln(  1.3  *  dtf  ) 

+  +  .1  *  sin<  1 .  +  2.5  *  dtf  ) 

+  +  .1  *  (  sln(  2.  +  3.5  *  dtf  ) 

+  +  sln(  3.  +  4.5  *  dtf  )  +  s 1 n (  4.  +  5.5  *  dtf  ) 

+  +  s  1  n  (  5.  +6.5*  dtf  )  +sln(  6.  +7.3*  dtf  )  )  ) 

g  3  anglo 
return 

end 

o 

subroutine  tgts(  J,  g  ) 

Implicit  double  precision  (  a— h,  o-z  ) 
common  /  tgts  /  a,  w 
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g  =  a  *  sln(  w  *  float!  j  )  ) 
return 
o 

end 

subroutine  fnmtn(n,x, fx, fun,e ,eps ,k) 

Implicit  double  precision  (  a-h,  o-z  ) 

dimension  xC 10) ,eps( 10) ,se(10),q(10),h(10,10),xl(10),xo(10) 
real#B  mj.lrada, 11 ,12,13,  lain,  mjfot 
o 

mjfot  »  2, 

o  reduced  from  20.  In  BRL  version  to  tame  subroutine 

m=n 

do  1  1  =  1  ,m 
se( 1 )=ops( 1 ) 
q(  1 ) =so ( 1 )#e 
xl ( 1 )=x( 1 ) 
xo ( 1 ) =xl ( 1 ) 
do  2  J=  1  ,m 

2  h(l , J)=0.0 
1  h( 1  , 1  )  =  1 .0 

o  1c  la  the  Iteration  counter  and  jc  is  the 

c  function  evaluation  counter. 

lc=  1 
Jc  =  0 
lr3=5 
go  to  112 

3  lmax=20*m 
fml n=  f bar 
f 0=  f bar 

f  J=  f bar 
del =0 .0 

assign  30  to  irl 

c  begin  Iteration 

50  do  41  J=  1  ,m 
ftJs<l(J) 

mj=  mjfct  *qj 
go  to  100 

30  q( J ) =max(so ( J ) , aba ( lmda )  ) 

1 f (abs (del ) .gt . abs ( f J-f bnr ) )go t o  41 
de 1  =  f  J-f bar 
Jd  =  J 

41  f  J=  f bar 

c  check  convergence 

l f ( 1 c .go . 1 max) go  to  91 
1  r2=  1 
kl  =  1 

psl 1=0.0 
emin=200. 
do  63  1  =  1  ,m 
t2=nbs(xl (1 )-xo( 1 ) ) 
lf(t2.oq.  0.)  go  to  63 
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if(t2,ge.eps(l)}  lr2=2 
pall =psl l+t2*t2 
t3=epa( 1 )/t2 
i f ( t3 . 1 t .emin)  emin=t3 
63  continue 

go  to  (90,70),ir2 

o  oheok  desirability  of  new  direction 

70  do  73  1  =  1, bj 

73  x( 1  )=xl ( 1 )+xl ( 1 )-xo( 1 ) 
ir3=6 
go  to  112 
75  fl»f0 

pat l=aqrt (pall) 
erainaemln*pai 1 
1 l=-pal 1 
f2=fmln 
12=0.0 
f3=f bar 
13=psl 1 

If (f3.ge.fl)goto  72 

1 f ( ( f l-(f2+f2)+f3)*( f l-f2-del )**2 .ge . . 5*del*( f l-f3)**2)go to  72 
o  compute  new  direotlon  and  use  directions 

o  i=l ,2,3, . . . , Jd-1 , Jd+1 , . . . ,n, new 

JJ=m-l 

1 f ( Jd-m)81 ,83,81 

81  do  82  l=Jd,JJ 
se ( i ) =so ( i+1 ) 
q(l)=q(l+l) 
do  82  Jl=l,m 

82  h(l,Jl)=h(i+i ,J1) 

83  do  84  Jl=l,m 

84  h(m, J1 )=(xi ( J1 )-xo( J1 ) )/psi 1 
se(m)=omln 

q(m) =psi 1 
qj=psi  1 

mj  =  mjfot  *  psi 1 
j=m 

assign  72  to  irl 
go  to  400 

o  proparo  for  new  iteration 

72  do  71  1  =  1  ,«n 

71  xo( 1 )=xl ( i  ) 
f0=fmln 

f J=fmin 
de  1  =0 
ic=lc+l 

assign  30  to  Irl 
go  to  50 

o  prepare  to  return 

91  kl=2 

90  do  92  i = 1 ,m 

92  x(i)=xl(i) 
fx=fmln 

If (k)93,96,93 


35 


.  i  .  ,,iv  ^  ■>,■■  ^v 


Jan  2  10:41 


93  k=kl 
97  return 

96  If (kl-l)94,97,94 

94  write(6,95)  lmax,kl 

93  format  (24h  funmln — not  converged — ,13,15k  Iterations, k  =,12) 
■  top 

o  find  minimum  along  a  line  (Initial  steps) 

100  12=0 

f2=fmln 
lmda=qj 
lr3=  1 
go  to  110 

102  if (fbar.gt.f2)goto  103 

11  =  12 

f  l  =  f2 
12=lmda 
f2=fbar 
lmda=qj+qj 
lr3=2 
go  to  110 
105  13= lmdn 
f3=f bar 
go  to  400 

103  13= lmda 
f3=fbar 
lmda=-qj 
1  r3=3 

go  to  110 

104  ll=lmda 
f l=fbar 

o  find  minimum  along  a  line 

400  t 1 = 12-13 
t2= 13-1 1 
13=11-12 
t4=  t I*t2*t3 
t5=  t l*f l+t2#f2+t3*f3 
t4=  t5/t4 
tl=ll*ll 
t2=12*12 
1 3= 1 3*13 

Imda=.5*((t2-t3)*fl+<t3-tl)*f2+<tl-t2)*f3)/t3 
If (t4)401 ,402,402 
401  If (abs(lmda)-mJ)403,4O3,4O2 

402  1 f ( f 1 . 1 t . f3)goto  404 
lmda=mj 

go  to  403 
404  lmda=-mj 

403  1 f (f 1 . 1 t . f2)goto  405 
If ( f 3 . 1 t .f2)goto  406 
lmln= 12 

fmln=f2 

407  1 f (abs( Imda-lml n) . 1 t . se ( J ) )  go  to  471 
If (lmda.eq.O.O)goto  408 

if (abs(( lmda- 1ml n)/ lmda) . 1 1 , .03) go  to  471 
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408  lr3-4 

go  to  110 

405  lmlnall 
fmln°f 1 
go  to  407 

406  lmlna13 
fralnaf3 
go  to  407 

480  if ( lmda.gt . 12)goto  481 
If ( lmda. 1 t . 1 1 )goto  482 
1 f ( f bar . 1 1 . f2)go to  483 

486  1 l3lmda 
f lafbar 
go  to  400 

481  1 f ( lmda.gt . 13)goto  484 
If (f bar. 1 t .f2)goto  483 

487  13a lmda 
f3=f bar 
go  to  400 

482  13=12 
f3=f2 
12=11 
f2=f  1 

go  to  486 

483  13=12 
f3=f2 

488  12= lmda 
f2=f bar 
go  to  400 

484  11=12 
fl  =  f2 
12=13 
f2=f3 

go  to  487 

485  11=12 
f  l  =  f2 

go  to  488 

471  linda=lmln 
fbar=fmln 
do  473  1 = 1 ,ra 

473  xi(i)=xl(i)+lmdn*li(J»l> 
go  to  irl.  (30,72) 

o  prepare  to  evaluate  function 

110  do  111  1=1, m 

111  x( 1 ) =xl ( 1 )+lmdn*h( J , 1 ) 

112  Jc=Jc+l 

f bar=fun(x,ro) 

c  special  for  control  ccccccco 
if  (  fbar  .le.  0.  )  go  to  91 
o  special  for  control  cccccccc 
go  to  (102,105,104,480,3,75  ),ir3 

end 
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1 


2 

3 


4 


c 


c 


subroutine  fft(  fr*  fi«  b  ) 


Implicit  ilfuble  preolsion 
dimension  fr(  1024  )»  fit 


(  a-ht  o-T. 
1024  ) 


) 


n  =  2  ##  k 


mr  =  0 
nn  3  n  -  1 


do  2  m  3  1 .  nn 
1  3  n 
1  =  1/2 

If  (  mr  +  1  .gft.  nn  )  go  to  l 
mr  =  mod(  rer*  1  )  +  1 
If  (  mr  .le.  m  )  go  to  2 
tr  3  fr(  m  +  1  ) 
fr(  m  +  1  )  3  frt  mr  +  1 
frt  mr  +  1  >  a  tr 
tl  3  fl<  m  +  1  ) 
fl (  m  +  1  )  3  f 1 (mr  +  l  > 
fit  mr  +  1  >  s  tl 


) 


continue 
1  3  1 

If  (  1  .ro.  n  )  return 
lutop  =  2  *  1 
el  =  1 

do  4  m  3  1 »  1  .  .  .  .  .  . 

n  =  3.1413926535  *  float!  1  -  m  )  /  el 

wr  3  cos (  a  ) 

vrl  3  sin!  a  ) 

do  4  t  3  m,  n.  Istep 
1=1  +  1 

t r  3  wr  #  fr(  J  1  -  wl  *  fl!  J  1 

tl  3  wr  *  fl!  J  >  +  "1  *  fr(  J  > 

frl  J  )  3  fr!  1  )  “  tr 

fl!  J  )  3  fl!  1  >  -  tl 

fr(l  )  53  frt  1  )  +  tr 

fl(l  >  s  f 1 1  1  )  +  1 1 

oont Inuo 
1  =  Istep 
ro  to  3 


end 

subroutine  manO!  delta,  nt ,  x  ) 


Implicit  double  precision  t  a-h,  o-x  > 

real*0  kk  „  , 

common  /  man  /  f?l  t  320  ),  r2,  rJ.  kk,  expt  , 

dimension  xt  10  ) 


tl  tl  , 


1 1  =  1 .  /  x !  I  ) 
kk  3  xt  2  ) 

1 1  =  xt  3  ) 
ttl  3  delta  /  tl 


1 1  tme 
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expt  =  exp(  —  1 1 1  ) 
tlti  =  tl  /  tl 
tltme  =  1.  -  tlti  -  expt 
do  50  n  =  1 ,  nt 
50  gl(  n  )  =  0. 

g2  =  0. 

g3  =  0. 
call  tgtaO 
return 
o 

end 

o 

subroutine  manl (  h,  J,  nt ,  diffdt,  dlfft,  gra  ) 
o 

Implicit  double  preolslon  (  a-h,  o-z  ) 
real*8  kk 

common  /  man  /  gl (  320  ),  g2,  g3,  kk,  oxpt ,  tlti,  tltme 
o 

do  60  n  =  1  ,  nt-1 
60  gl  (  n  )  =  gl  (  n+1  ) 

call  tgtal(  J,  dither,  g  > 
dlfft  =  g  -  h 
diffdt  =  g  +  dither  -  h 
gl (  nt  )  =  diffdt 

g21  =  g2 

B 2  =  kk  *  gl (  1  ) 

g3  =  oxpt  *  g3  +  tltrao  *  g21  +  tlti  *  g2 
gm  =  g3 
return 
c 

end 

subroutine  plantOt  b,  kb,  bk,  bkl ,  cf,  ds,  dt  ) 

o 

Implicit  double  precision  (  a~h,  o-z  ) 
re a 1*8  JJ,  kb 

common  /  plant  /  expb,  bdxpb,  bklsh,  bklshl,  efdt,  dedspc, 
+  delta,  g3b,  g4,  hn ,  hnbl 

o 

JJ  =  1. 

delta  =  dt 
t J  =  delta/  JJ 
expb  =  exp(  -b  *  kb  *  tj  ) 
bdxpb  =  (  1 .  -  expb  )  /  b 
bklsh  =  bk 
bklshl  =  bkl 
cfdt  =  cf  *  tj 
dedspc  =  ds 
return 
o 

end 

o 

subroutine  plant 1(  h  ) 
c 
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Implicit  double  precision  (  a-h,  o-z  ) 

common  /  plant  /  expb,  bdxpb,  bklsh,  bklahl ,  ofdt,  dedapo, 
+  delta,  g3b,  g4,  hn,  hnbl 
o 

g3b  =  0. 
g4  =  0. 
h  =  0. 
hn  =  0 . 
hnbl  =  0. 
return 
o 

end 

c 

subroutine  plant2(  g3,  h  > 
c 

Implicit  double  precision  (  a-h,  o-z  ) 

common  /  plant  /  expb,  bdxpb,  bklsh,  bklahl,  ofdt,  dedspo, 
+  delta,  g3b,  g4 ,  hn,  hnbl 
o 

g3bl  =  g3b 
g3b  =  0. 

If  (  abs(  g3  )  ,gt.  dedspo  > 

+  g3b  =  g3  -  slgn(  dedspc ,  g3  ) 
c  If  (  abs(  g3b  -  g3bl  )  ,gt.  bklshl  ) 

o  +  g3b  =  g3b  -  slgn(  bklshl,  g3b  -  g3bl  ) 
g4  =  expb  *  g4  +  bdxpb  #  g3bl 
g4  =  g4  -  slgn(  mini  ofdt,  aba(  g4  )  ) ,  g4  ) 
hnbl  =  hnbl  +  delta  *  g4 
If  (  absl  hnbl  -  hn  )  ,gt.  bklsh  ) 

+  hn  =  hnbl  -  sign(  bklsh,  hnbl  -  hn  ) 
h  =  hn 
re  turn 
o 

end 

subroutine  statsOt  dt ,  Jmax  ) 

c 

Implicit  double  preolslon  (  a-h,  o-z  ) 
common  /  stats  /  delta,  points,  tmax, 

+  aura,  sumsq,  sumin,  suraout ,  suml2,  sumo2,  < 

+  sunt ,  sumt2,  sumlt,  sumot 
c 

delta  -  dt 

points  =  jmax 

tmax  =  points  *  delta 

sum  =  0 . 

sumsq  =  0. 

sumin  =  0. 

sumout  =  0. 

sum 12  =  G. 

sum02  =  0. 

sumt  =  0. 

surat 2  =  0. 

sumlt  =0. 
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sumo  t  =  0, 
return 
c 

end 

c 

subroutine  statsl  (  dlfft,  frln,  frout,  J  ) 
c 

implicit  double  precision  (  a-h,  o-z  ) 
common  /  stats  /  delta,  points,  tmax, 

+  sum,  sumsq,  sumin,  sumout,  suml2,  sumo2, 

+  sunt ,  sumt2,  sumlt,  sumot 
o 

dt  =  dlfft 

sum  =  sum  +  dt  *  delta 

sumsq  =  sumsq  +  dt  *  dt  *  delta 

sumin  =  sumin  +  frin 

sumout  =  sumout  +  frout 

sumi2  =  sum!2  +  frln  *  frin 

sumo2  =  sumo2  +  frout  *  frout 

tj  =  J 

sumt  =  sumt  +  tj 
sumt2  =  sumt2  +  tj  *  tj 
sumlt  =  sumlt  +  frln  *  tj 
sumot  =  sumot  +  frout  *  tj 
return 
o 

end 

o 

stibrout  Ine  stats2  (  nvln,  avout ,  rmsd  ) 
c 

Implicit  double  precision  (  a-h,  o-z  ) 
common  /  stats  /  delta,  points,  troax, 

+  sum,  sumsq,  sumin,  sumout,  sural2,  sumo2, 

+  sumt,  sumt2,  sumlt,  sumot 

common  /  statw  /  rms ,  aveln,  aveout , 

+  devln,  devout,  devid,  dovod,  ai ,  bl ,  ao,  bo 
o 

rms  =  sqrt(  sumsq  /  tmax  ) 
rmsd  =  rms 

aveln  =  sumin  /  points 

uv  1  n  -  ave  1  n 

aveout  =  sumout  /  points 

avout  =  aveout 

sl2  =  sum!2  -  sumin  %  aveln 

so2  =  sumo2  -  sumout  *  aveout 

ptsl  =  points  -  1 . 

devln  =  sqrt(  s!2  /  ptsl  ) 

devout  =  eqrt(  so2  /ptsl  ) 

avot  =  sumt  /  points 

st2  =  sumt 2  -  sumt  *  avet 

sit  =  sumlt  -  sumin  *  avot 

sot  =  sumot  -  sumout  *  avet 

bl  =  sit  /  st2 

bo  =  sot  /  st2 
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al  =  aveln  -  bl  *  avet 
ao  =  aveout  -  bo  *  avet 
pts2  =  points  -  2 

devid  =  sqrt  (  (  si2  -  bl  *  si t  >  /  pts2  ) 
dovod  =  sqrt  (  (  so2  -  to  *  sot  )  /  pts2  ) 
return 

end 

subroutine  statsw 

implicit  double  preeision  (  a-h,  o-z  ) 


common 

/  statw 

/ 

rras  ♦ 

aveln,  aveout, 

devln, 

devout 

+  devid, 

devod, 

ai 

,  bl 

,  ao,  bo 

write  ( 

6,  100 

) 

rras , 

aveln,  aveout , 

devln, 

devout 

+  devid, 

devod , 

al 

♦  bi 

,  ao,  bo 

100  format  (  1  RMS  AVEIN  AVEOUT  DEVIN  DEVOUT  ’ 
+  ’DEVID  DEVOD  AI  BI  AO  BO  ’ 

+  /  3pllf6.3  /  ) 

return 
c 

end 
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